%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This code is written to estimate the GARCH-MIDAS model with rolling window RV
%
% This code needs the data of the following format:
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% stkmkt_ret_data
% yyyy  mm  dd  yyyymmdd  vwretd  ewretd
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% For practical purposes, the parameter m and theta in the tau function
% are squared as follows:
%
% tau = (m^2) + (theta^2) * ( X * betapolyn(K,[(K-1):-1:1]',w) )
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


clear all; 
close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% daily stock returns: 1926 - 2007 / NYSE-AMEX-NASDAQ
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% stkmkt_ret_data
% yyyy  mm  dd  yyyymmdd  vwretd  ewretd
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp('which spanning horizon?')
term=input('(monthly=1/quarterly=2/biannual=3/annual=4)    ');

disp('   ')
K=input('how many daily lags for MIDAS filter?    ');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We are using daily returns of stocks listed in NYSE/AMEX/NASDAQ
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load stkmkt_ret_data;
daily_stk = [stkmkt_ret_data(:,1:3) stkmkt_ret_data(:,5)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[T0,C0]=size(daily_stk);

ft_ind0=find( diff(daily_stk(:,2)) ~= 0 ); % row number
ft_ind0=ft_ind0+1; % row number where new month begins

ft_ind=zeros(T0,1);
ft_ind(ft_ind0,1)=1; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
daily_stk=[daily_stk(:,1:3) ft_ind daily_stk(:,4)];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ind1=find( diff(daily_stk(:,2)) ~= 0 ); % row number
ind1=ind1+1; %

num_days = zeros(T0,1);
num_days(ind1)=[diff(ind1);(T0+1)-ind1(end)];

num_days(1) = ind1(1)-1;

%--------------------------------------------------------------
daily_stk=[daily_stk(:,1:4) num_days daily_stk(:,5)];
%--------------------------------------------------------------

        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% realized volatility (rolling)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
switch term
    case 1
        spanning=22;
    case 2
        spanning=65;
    case 3
        spanning=125;
    case 4
        spanning=250;
end;
        
rolling_RV=zeros(T0,1);
for i=(spanning+1):T0
     rolling_RV(i)=sum( daily_stk(i-spanning:i-1,6).^2 );
end;
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
daily_stk = [daily_stk rolling_RV];
daily_stk = daily_stk(spanning+1:end,:); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The first spanning number of days of data are burned.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (1) yr 
% (2) mth 
% (3) day 
% (4) 'beginning of new month index (0 or 1)' 
% (5) num_days in month
% (6) ret
% (7) rolling RV


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subsample choice
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

date=daily_stk(:,1)*10000+daily_stk(:,2)*100+daily_stk(:,3);

disp('    ')
disp('Sample Period? (ex. [19270101;20071231] / for whole sample, just press enter)');
disp('[The first K days (the number of the daily lags in MIDAS)') 
disp('of data from your first date will be burned]')
subperiod=input('Your Choice?  ');

if isempty(subperiod)==1
    subperiod=[date(1);date(end)];
else
end;

start_ind=find(date==subperiod(1));
j=1;
while isempty(start_ind)==1
    start_ind=find(date==(subperiod(1)+j));
    j=j+1;
end;

end_ind=find(date==subperiod(2));
j=1;
while isempty(end_ind)==1
    end_ind=find(date==(subperiod(2)-j));
    j=j+1;
end;

daily_stk=daily_stk(start_ind:end_ind,:);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r = daily_stk(K:end,6);   

disp('starting date for fitted return series data')
daily_stk(K,1:3)

start_yr  = daily_stk(K,1);
end_yr    = daily_stk(end,1);


% Constructing Quarterly RV series -------------------------

qtr_ind_temp = find( ((daily_stk(K:end,2)==1)|(daily_stk(K:end,2)==4)|(daily_stk(K:end,2)==7)|(daily_stk(K:end,2)==10))&(daily_stk(K:end,4)==1) );
if qtr_ind_temp(1)~=1
    qtr_ind_temp = [1;qtr_ind_temp];
else
end;
n_qtr_ind_temp =length(qtr_ind_temp);
qtr_ind = (K-1) + qtr_ind_temp;

N_qtr  = length(qtr_ind);
RV_qtr = zeros(N_qtr,1);
Quarticity = zeros(N_qtr,1);

for i=1:(N_qtr-1)
    RV_qtr(i) = sum( daily_stk( qtr_ind(i):(qtr_ind(i+1)-1) , 6 ).^2 );
    Quarticity(i) = sum( daily_stk( qtr_ind(i):(qtr_ind(i+1)-1) , 6 ).^4 );
end;

RV_qtr(end) = sum( daily_stk( qtr_ind(N_qtr):end , 6 ).^2 );
%--------------------------------------------------------


T=length(r);
X=zeros(T,K-1);

for i=1:T
    X(i,:)=daily_stk((i+1):K+(i-1),7)';   
end;

r_yr = daily_stk(K:end,1);

ref_yr = [start_yr+(5-mod(start_yr,5)):5:end_yr - mod(end_yr,5)]';

if mod(start_yr,5)==0
    ref_yr=[start_yr;ref_yr];
else
end;

ref_yr_length=length(ref_yr);
tick_yr=zeros(ref_yr_length,1);

for i=1:ref_yr_length
    tick_yr(i)=min(find(r_yr==ref_yr(i)));
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ESTIMATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% The following are the starting values: 

miu_i    =0.0001;
alpha_i  =0.10;
beta_i   =0.85;
theta_i  =0.1;
w_i      =10;
m_i      =0.001;

parameter_i=[miu_i;alpha_i;beta_i;theta_i;w_i;m_i];

g1=1;

options=optimset('Display','iter','MaxFunEvals',3000,'MaxIter',3000,'TolFun',1e-6,'TolX',1e-6);
options =optimset(options,'GradObj','off');
options=optimset(options,'DerivativeCheck','off');
options=optimset(options,'Diagnostics','off');

A=[
    0  1  1  0  0  0;
   -1  0  0  0  0  0;
    0 -1  0  0  0  0;
    0  0 -1  0  0  0;
    0  0  0 -1  0  0;
    0  0  0  0 -1  0;
    0  0  0  0  1  0;
    0  0  0  0  0 -1
  ];

b=[
      .99999;
           0;
           0; 
           0;
           0;
           0;
         100;
  -0.0000001;
   ];


disp('Optimization Results ------------------')
[parameter,fval,exitflag,output,lambda,grad,hessian] = fmincon('nlh_garchmidas',parameter_i,A,b,[],[],[],[],[],options,g1,K,r,X);
disp('---------------------------------------')

T_n=length(r);

Cov_Var=hessian/T_n;
V=inv(Cov_Var);
S_var=diag(V)/T_n;
S_err=sqrt(S_var);

LLF = -nlh_garchmidas(parameter,g1,K,r,X);
BIC = -2*LLF/T_n + 6*log(T_n)/T_n;


disp(['Likelihood= ' num2str(-nlh_garchmidas(parameter,g1,K,r,X))])
disp(' ')

name=['miu=   ';'alpha= ';'beta=  ';'theta= ';'w=     ';'m=     '];

disp('name      estimates        ')
disp([name num2str(parameter)])
disp('   ')
disp('name      s.e.        ')
disp([name num2str(S_err)])
disp('   ')
disp('name      t-stat        ')
tstat=parameter./S_err;
disp([name num2str(tstat)])


miu   =parameter(1);
alpha =parameter(2);
beta  =parameter(3);
theta =parameter(4);
w     =parameter(5);
m     =parameter(6);


tau = (m^2) + (theta^2) * ( X * betapolyn(K,[(K-1):-1:1]',w) ); % function of theta and w 

g=zeros(T_n,1);
g(1)=g1;

for i=2:T_n
    g(i)=(1-alpha-beta)+alpha*((r(i-1)-miu)^2)/tau(i-1)+beta*g(i-1);
end;


e=(r-miu)./sqrt(tau.*g);

figure
hist(e,100);title('Distribution of errors')

figure
plot(flipdim(betapolyn(K,[(K-1):-1:1]',w),1));title('optimal MIDAS weights')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% var_qtr = zeros(N_qtr,1);
% tau_qtr = zeros(N_qtr,1);
% 
% for i=1:N_qtr-1
%     var_qtr(i) = sum( tau( (qtr_ind(i)-(K-1)) : (qtr_ind(i+1)-1-(K-1)) ) .* g((qtr_ind(i)-(K-1)):(qtr_ind(i+1)-1-(K-1))) );
%     tau_qtr(i) = sum( tau( (qtr_ind(i)-(K-1)) : (qtr_ind(i+1)-1-(K-1)) ) );
% end;
% 
% var_qtr(end) = sum( tau((qtr_ind(i)-(K-1)):(qtr_ind(i+1)-1-(K-1))) .* g((qtr_ind(i)-(K-1)):(qtr_ind(i+1)-1-(K-1))) );
% tau_qtr(end) = sum( tau((qtr_ind(i)-(K-1)):(qtr_ind(i+1)-1-(K-1))) );
% 
% qtrRV_series  = zeros(T_n,1);
% qtrVar_series = zeros(T_n,1);
% 
% qtr_I = zeros(T_n,1);
% qtr_I(qtr_ind-(K-1)) = 1;
% p = 0;
% 
% for i=1:T_n
%     p = p + qtr_I(i);
%     qtrRV_series(i)  = RV_qtr(p);
%     qtrVar_series(i) = var_qtr(p);
% end;


figure
% subplot(2,1,1)
plot(sqrt(252*tau.*g),'g--','LineWidth',1);
hold on
plot(sqrt(252*tau),'b-','LineWidth',2);
axis([0 T_n 0.05 0.5])
legend('(ann.) conditional volatility (\tau*g)^{1/2}','(ann.) secular component (\tau)^{1/2}')
title('conditional volatility and its long run component of stock market returns');
xlabel('year')
ylabel('annualized volatility')
set(gca,'XTick',tick_yr);
set(gca,'XTickLabel',num2str(ref_yr));

% subplot(2,1,2)
% plot(sqrt(4*RV_qtr),'Color',[0.6 0.6 0.6],'LineStyle',':','Marker','s','MarkerSize',3,'LineWidth',1)
% hold on
% plot(sqrt(4*var_qtr),'g--','LineWidth',2)
% plot(sqrt(4*tau_qtr),'b-','LineWidth',2)
% axis([0 N_qtr 0.02 0.5])
% legend('(ann.) quarterly RV^{1/2}','(ann.) quarterly aggregated \tau*g^{1/2}','(ann.) quaterly aggregated \tau^{1/2}')
% title('quarterly aggregation of \tau*g / \tau and quarterly RV');
% xlabel('year')
% ylabel('annualized volatility')
% 
% r_yr3 = [
%           repmat(daily_stk(K,1),(5-(floor((daily_stk(K,2)-1)/3)+1)),1);
%           reshape(repmat([(daily_stk(K,1)+1):1:daily_stk(end,1)],4,1),4*(daily_stk(end,1)-daily_stk(K,1)) ,1)
%         ];
% 
% 
% for i=1:ref_yr_length
%     tick_yr3(i)=min(find(r_yr3==ref_yr(i)));
% end;
% 
% 
% set(gca,'XTick',tick_yr3);
% set(gca,'XTickLabel',num2str(ref_yr));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

